in Nonmetropolitan than Metropolitan Area

An Overview of the Data

This interactive plot examines the differences of Standardized Mortality Ratio (SMR), the proportion between observed death and expected death, of the five causes of death in each public health region. Each observation indicates for which state it belongs to with the size of which is reflective of the corresponding population. Here we compare SMR within the same geographic area. The differences of the distribution of the SMR among five causes of deaths might associate with the sociodemographic, cultural, behavioral and structural factors of the specific region.

cod_data %>%
  mutate(year = as.factor(year))
## # A tibble: 191,748 x 12
##      year age_range cause_of_death   state        locality observed_deaths
##    <fctr>     <chr>          <chr>   <chr>           <chr>           <int>
##  1   2005      0-49         Cancer Alabama             All             756
##  2   2005      0-49         Cancer Alabama    Metropolitan             556
##  3   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  4   2005      0-49         Cancer Alabama             All             756
##  5   2005      0-49         Cancer Alabama    Metropolitan             556
##  6   2005      0-49         Cancer Alabama Nonmetropolitan             200
##  7   2005      0-49         Cancer Alabama             All             756
##  8   2005      0-49         Cancer Alabama    Metropolitan             556
##  9   2005      0-49         Cancer Alabama Nonmetropolitan             200
## 10   2005      0-54         Cancer Alabama             All            1346
## # ... with 191,738 more rows, and 6 more variables: population <int>,
## #   expected_deaths <int>, potentially_excess_deaths <int>,
## #   percent_excess_death <dbl>, mortality <dbl>, hhs_region <int>
p <- cod_data %>%
  plot_ly(
    x = ~expected_deaths, 
    y = ~observed_deaths, 
    size = ~population, 
    color = ~cause_of_death, 
    frame = ~hhs_region, 
    text = ~state, 
    hoverinfo = "text",
    type = 'scatter',
    mode = 'markers'
  ) %>% 
  layout(title = "Change of Standardized Mortality Ratio in National Public Health Regions")

p

Here we show the division of the 10 U.S. Department of Health and Human Services public health regions in the figure below.

From the linear regression plot below, we can see that the mortality of cancer, unintentional injury in both metropolitan regions and nonmetropolitan regions decrease through years 2005-2015. In contrast, the mortality of heart disease and chronic lower respiratory disease in both metropolitan regions and nonmetropolitan regions increase through years 2005-2015. For stroke, the mortality in metropolitan regions decreases through years 2005-2015 whereas the mortality in metropolitan regions slightly increases through years 2005-2015. Concerning mortality and the trend in mortality, unintentional injury and heart disease needs more attention from public health experts.

#mortality~year+cause
cod_datayr2 = read_csv("./data/NCHS_-_Potentially_Excess_Deaths_from_the_Five_Leading_Causes_of_Death.csv") %>%
  clean_names() %>%
  na.omit() %>%
  filter(!(state == "United States")) %>%
  separate(., percent_potentially_excess_deaths, into = c("percent_excess_death"), sep = "%") %>%
  mutate(percent_excess_death = as.numeric(percent_excess_death), mortality = observed_deaths/population * 10000, mortality = as.numeric(mortality)) %>%
  select(year, cause_of_death, state, benchmark, locality, observed_deaths, population, expected_deaths, potentially_excess_deaths, percent_excess_death, mortality)
## Parsed with column specification:
## cols(
##   Year = col_integer(),
##   `Cause of Death` = col_character(),
##   State = col_character(),
##   `State FIPS Code` = col_character(),
##   `HHS Region` = col_integer(),
##   `Age Range` = col_character(),
##   Benchmark = col_character(),
##   Locality = col_character(),
##   `Observed Deaths` = col_integer(),
##   Population = col_integer(),
##   `Expected Deaths` = col_integer(),
##   `Potentially Excess Deaths` = col_integer(),
##   `Percent Potentially Excess Deaths` = col_character()
## )
## Warning: Too many values at 191748 locations: 1, 2, 3, 4, 5, 6, 7, 8, 9,
## 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, ...

lineplotyr = cod_datayr2%>%
  filter(benchmark == "2010 Fixed",
         locality == c("Metropolitan", "Nonmetropolitan"))%>%
  group_by(cause_of_death,year, locality)%>%
  distinct(year, cause_of_death, locality, .keep_all = TRUE)%>%
  ungroup(year, cause_of_death)%>%
  ggplot(aes(x = year, y = mortality, group = locality, color = locality))+
  geom_point(size = 0.5)+
  geom_smooth(method = lm, size = 0.5,se=FALSE) +
    facet_grid(.~cause_of_death)+
  theme(axis.text.x = element_text(angle = 45, hjust = 1))+
  scale_x_continuous( breaks=seq(2005,2015,1))
ggplotly(lineplotyr) %>% 
  layout(title = "Mortality over years for five leading causes of death in different locality",
     #xaxis = list(showticklabels = FALSE),
     legend = list(orientation = "h",
                   y = -0.1, x = 0.3))
## We recommend that you use the dev version of ggplot2 with `ggplotly()`
## Install it with: `devtools::install_github('hadley/ggplot2')`

Descriptive Analysis

cod_data %>%
  mutate(cause_of_death = fct_reorder(cause_of_death, mortality)) %>% 
  group_by(cause_of_death) %>% 
  ggplot(aes(x = cause_of_death, y = mortality, color = cause_of_death)) +
  geom_boxplot() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title="Locality vs. Mortality") +
  labs(x="Locality", y="Mortality", color = "Cause of Death") +
  facet_grid(~locality)

This bar graph shows the distribution of cause of death within mortality in the three geographic regions: Metropolitan, Nonmetropolitan and All. We observe that the rank in causes of death is the same regardless of the locality. Nonmetropolitan area seems to have the highest mortaility, the number of deathes observed in every 10000 people. We observe that the median of mortaility of Cancer is the highest among all 5 causes of death, followed by heart disease, unintentional injury, heart disease, chronic lowee respiratory disease and stroke. All 5 causes of death are right skewed regardless of the locality.

Linear Model

To further support our graphs above, linear regression is performed between Mortality and Locality. Here we obtain the linear regression model of mortality = 4.103 - 0.226localityMetropolitan -1.145localityNonmetropolitan. We expect to see 0.22576 less deaths for every 10000 people in the metropolitan region as compared to all regions. We expect to see 1.14451 more death for every 10000 people in the nonmetropolitan region as compared to all regions.

mortality_lm = lm(mortality ~locality, data = cod_data)
summary(mortality_lm)
## 
## Call:
## lm(formula = mortality ~ locality, data = cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.1454 -2.9547 -1.2983  0.9988 19.8959 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              4.10331    0.01672 245.474   <2e-16 ***
## localityMetropolitan    -0.22576    0.02372  -9.519   <2e-16 ***
## localityNonmetropolitan  1.14451    0.02438  46.942   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.317 on 191745 degrees of freedom
## Multiple R-squared:  0.01821,    Adjusted R-squared:  0.0182 
## F-statistic:  1778 on 2 and 191745 DF,  p-value: < 2.2e-16
rl3 = broom::tidy(mortality_lm)
kable(rl3)

Data Visualization

Does mean percent excess death have significant difference between metro and non-metro groups?

gp_cod_data = cod_data %>%
  group_by(year, locality) %>% 
  summarise(mean_ped = mean(percent_excess_death)) 

gp_cod_data_Metro = gp_cod_data %>%
  group_by(year) %>%
  filter(locality == "Metropolitan")
gp_cod_data_nonMetro = gp_cod_data %>%
  group_by(year) %>%
  filter(locality == "Nonmetropolitan")
gp_cod = cbind(gp_cod_data_Metro, gp_cod_data_nonMetro)

p = gp_cod %>%
plot_ly( x = ~year, y = ~mean_ped, name = 'Metropolitan', type = 'scatter', mode = 'markers') %>%
  add_trace( y = ~mean_ped1, name = 'Nonmetropolitan', mode = 'markers') %>% 
  layout(yaxis = list(title = 'mean percent excess death (Metropolitan/Nonmetropolitan)'))
p

The above plot shows the gap of mean percent excess death for metropolitan and nonmetropolitan. Nonmetropolitan has a obviously larger mean percent excess death than metropolitan areas. In order to see if the difference is significant, we will run a regression analysis regarding mean excess death and locality.

Regression Analysis

The regression model is mean percent_excess_death = 33.89 - 2.05Metropolitan + 8.046Nonmetropolitan.
For people in metropolitan, the expected percent excess death has a decrease of 2.06. * For people in non-metropolitan aream the expected percent excess death has an increase of 8.046. * The p-values show these predictors are both significant. The adjusted R squared is 0.91, which indicates the regression model explains 91% of variation in mean percent excess death due to the variation in locality.

gp_locality_lm = lm(mean_ped ~locality , data = gp_cod_data)
summary(gp_locality_lm)
## 
## Call:
## lm(formula = mean_ped ~ locality, data = gp_cod_data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.5667 -1.0325 -0.3470  0.7922  3.2937 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)              33.8866     0.4071  83.238  < 2e-16 ***
## localityMetropolitan     -2.0599     0.5757  -3.578   0.0012 ** 
## localityNonmetropolitan   8.0460     0.5757  13.975 1.13e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.35 on 30 degrees of freedom
## Multiple R-squared:  0.9198, Adjusted R-squared:  0.9145 
## F-statistic: 172.1 on 2 and 30 DF,  p-value: < 2.2e-16
rl4 = broom::tidy(gp_locality_lm)
kable(rl4)

The regression model is mean percent_excess_death = 33.89 - 2.05*Metropolitan + 8.046Nonmetropolitan. For people in metropolitan, the expected percent excess death has a decrease of 2.06. For people in non-metropolitan aream the expected percent excess death has an increase of 8.046. The p-values show these predictors are both significant. The adjusted R squared is 0.91, which indicates the regression model explains 91% of variation in mean percent excess death due to the variation in locality.

=======

3d938490d6122be3e5a56fec6840c636cf3df2be

Account for Locality

For each of the five cause of death, what is the mean percent of excess death? It is different across locality.

As shown in the boxplots for all three regions, the rank of Mean Percentage of Excess Death of the five cause of death remain the same for different localities. For all five causes, the mean percent of death is lower for Metropolitan and higher for Nonmetropolotan.

cod_data %>%
  select(state, locality, percent_excess_death, hhs_region, cause_of_death) %>%
  filter(locality != "All") %>% 
  mutate(hhs_region = as.factor(hhs_region)) %>% 
  group_by(cause_of_death, locality) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  ungroup() %>% 
  ggplot(aes(x = cause_of_death,  y = mean_ped, fill = cause_of_death)) +
  geom_col() +
  facet_grid(~locality) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "Five Leading Causes of Death",
    y = "Mean Percentage of Excess Death",
    fill = "Causes of Death"
  )

Account for HHS Regions

region_lm = lm(region_cod_data$mean_ped~region_cod_data$hhs_region)
rl1 = broom::tidy(region_lm)
kable(rl1)

region_locality_lm = lm(region_cod_data$mean_ped ~ region_cod_data$hhs_region+region_cod_data$locality)
rl2 = broom::tidy(region_locality_lm)
kable(rl2)

cod_data %>%
  select(state, locality, percent_excess_death, hhs_region, cause_of_death) %>%
  filter(locality != "All") %>% 
  group_by(cause_of_death, locality, hhs_region) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  mutate(hhs_region = as.factor(hhs_region)) %>% 
  ggplot(aes(x = hhs_region,  y = mean_ped, fill = hhs_region)) +
  geom_col() +
  facet_grid(~cause_of_death) +
  theme(axis.text.x = element_text(angle = 15)) +
    labs(
      y = "Mean %Excess Death",
      x = "HHS Regions",
      fill = "HHS Region"
    ) 

map_cod_data = cod_data %>%
  filter(locality == "Metropolitan") %>%
  select(state, locality, percent_excess_death) %>% 
  group_by(state) %>% 
  summarise(mean_ped = mean(percent_excess_death)) %>% 
  dplyr::filter(!(state == "District of\nColumbia"))
  

map = as.tibble(fifty_states) %>%
  group_by(id) %>% 
  summarize(clong = mean(long), clat = mean(lat)) %>% 
  filter(!(id == "district of columbia"))

df <- cbind(map, state.abb, state.center, rate = unique(map_cod_data$mean_ped))

ggplot(df, aes(map_id = id)) + 
    geom_map(aes(fill = rate), map = fifty_states) + 
    expand_limits(x = fifty_states$long, y = fifty_states$lat) + 
    labs(x = "", y = "") +
    theme(panel.background = element_blank(), 
          axis.text.x = element_blank(), 
          axis.text.y = element_blank(), 
          axis.ticks = element_blank()) + 
    geom_text(aes(x = clong, y = clat, label = state.abb)) +
  scale_fill_gradient(low="gold", high="red")